function matr = matrix_of_boundary_conditions(p, m, CONSTS, varargin)

    if (size(varargin, 1) >= 1)
        use_combination = varargin{1, 1};
    else
        use_combination = false;
    end

    k0    = CONSTS.k0;
    a     = CONSTS.a;
    eta   = CONSTS.eta;
    eps_a = CONSTS.eps_a;

    s  = func_s(p, CONSTS);
    q1 = func_q1_from_p(p, CONSTS);
    q2 = func_q2_from_p(p, CONSTS);
    Q1 = k0*a*q1;
    Q2 = k0*a*q2;
    S  = k0*a*s;
    
    bessel_m_q1 = besselj(m, Q1);
    bessel_m_q2 = besselj(m, Q2);
    bessel_K_m = besselk(m, S);

    n1     = func_n1_from_p(p, CONSTS);
    n2     = func_n2_from_p(p, CONSTS);
    K_m    = func_K_m(p, m, CONSTS);
    J_m1   = func_J_m1(p, m, CONSTS);
    J_m2   = func_J_m2(p, m, CONSTS);
    J_m1_1 = func_J_m1_1(p, m, CONSTS);
    J_m2_1 = func_J_m2_1(p, m, CONSTS);
    
    if (use_combination)
        matr = [1.0,         1.0,       0.0,                              -1.0;  ...
                n1,          n2,        -(eta/eps_a),                     0;     ...
                J_m1,        J_m2_1,    -(p./eps_a).*(m./((k0*a*s).^2)),  K_m;   ...
                n1.*J_m1_1,  n2.*J_m2,  K_m,                              -(p.*(m./((k0*a*s).^2)))];
    else
        matr = [Q1.*bessel_m_q1,              Q2.*bessel_m_q2,            0.0,                                     S.*bessel_K_m;  ...
                n1.*Q1.*bessel_m_q1,          n2.*Q2.*bessel_m_q2,        (eta/eps_a)*S.*bessel_K_m,               0;     ...
                J_m1.*Q1.*bessel_m_q1,        J_m2_1.*Q2.*bessel_m_q2,    (p./eps_a).*(m./(S.^2)).*S.*bessel_K_m,  -K_m.*S.*bessel_K_m;   ...
                n1.*J_m1_1.*Q1.*bessel_m_q1,  n2.*J_m2.*Q2.*bessel_m_q2,  -K_m.*S.*bessel_K_m,                     (p.*(m./(S.^2))).*S.*bessel_K_m];
    end
end

function y = func_J_m1(p, m, CONSTS)
    k0 = CONSTS.k0;
    a = CONSTS.a;
    q1 = func_q1_from_p(p, CONSTS);
    alfa_1 = func_alfa_1(p, CONSTS);
    bessel_m_q1 = besselmx(double('J'), m, k0*a*q1, 0);
    bessel_mp1_q1 = besselmx(double('J'), m+1, k0*a*q1, 0);

    y = (bessel_mp1_q1./((k0*a*q1).*bessel_m_q1)) + ...
           ((alfa_1.*m)./((k0*a*q1).^2));
end

function y = func_J_m2(p, m, CONSTS)
    k0 = CONSTS.k0;
    a = CONSTS.a;
    q2 = func_q2_from_p(p, CONSTS);
    beta_2 = func_beta_2(p, CONSTS);
    bessel_m_q2 = besselmx(double('J'), m, k0*a*q2, 0);
    bessel_mp1_q2 = besselmx(double('J'), m+1, k0*a*q2, 0);

    y = (bessel_mp1_q2./((k0*a*q2).*bessel_m_q2)) - ...
           ((beta_2.*m)./((k0*a*q2).^2));
end

function y = func_J_m1_1(p, m, CONSTS)
    k0 = CONSTS.k0;
    a = CONSTS.a;
    q1 = func_q1_from_p(p, CONSTS);
    beta_1 = func_beta_1(p, CONSTS);
    bessel_m_q1 = besselmx(double('J'), m, k0*a*q1, 0);
    bessel_mp1_q1 = besselmx(double('J'), m+1, k0*a*q1, 0);

    y = (bessel_mp1_q1./((k0*a*q1).*bessel_m_q1)) - ...
           ((beta_1.*m)./((k0*a*q1).^2));
end

function y = func_J_m2_1(p, m, CONSTS)
    k0 = CONSTS.k0;
    a = CONSTS.a;
    q2 = func_q2_from_p(p, CONSTS);
    alfa_2 = func_alfa_2(p, CONSTS);
    bessel_m_q2 = besselmx(double('J'), m, k0*a*q2, 0);
    bessel_mp1_q2 = besselmx(double('J'), m+1, k0*a*q2, 0);

    y = (bessel_mp1_q2./((k0*a*q2).*bessel_m_q2)) + ...
           ((alfa_2.*m)./((k0*a*q2).^2));
end

function y = func_K_m(p, m, CONSTS)
    k0 = CONSTS.k0;
    a = CONSTS.a;
    s = func_s(p, CONSTS);

    y = (besselk(m+1, k0*a*s, 0)./((k0*a*s).*besselk(m, k0*a*s, 0))) - ...
           (m./((k0*a*s).^2));
end
